1 Wstęp

Będziemy badać metodę svm (zaimplementowaną w pakiecie e1071) na zbiorach apartments z pakietu DALEX, na którym będziemy dokonywać regresji, oraz na zbiorze bodyfat z bazy zbiorów OpenML, na której będziemy dokonywać klasyfikacji.

bodyfat <- getOMLDataSet(data.name = "bodyfat")
bodyfat <- bodyfat$data

len1 <- nrow(apartments)
len2 <- nrow(bodyfat)

trainind1 <- sample(len1, ceiling(len1*6/10))
testind1 <- setdiff(1:len1, trainind1)
len1tr <- length(trainind1)
len1ts <- length(testind1)

trainind2 <- sample(len2, ceiling(len2*6/10))
testind2 <- setdiff(1:len2, trainind2)
len2tr <- length(trainind2)
len2ts <- length(testind2)

2 Dopasowanie modeli

Na początek po prostu dopasujemy modele do danych. Zrobimy to z wykorzystaniem frameworka mlr.

task1 <- makeRegrTask("regr_apartments", data = apartments[trainind1,], target = "m2.price")
task2 <- makeClassifTask("classif_cosinnegos", data = bodyfat[trainind2,], target = "binaryClass")

lrn1 <- makeLearner("regr.svm", par.vals = list(scale = FALSE))
lrn2 <- makeLearner("classif.svm", predict.type = "prob", par.vals = list(scale = FALSE))

model1 <- train(lrn1, task1)
model2 <- train(lrn2, task2)

out1 <- predict(model1, newdata = apartments[testind1,])
out2 <- predict(model2, newdata = bodyfat[testind2,])

Wyniki regresji na apartments:

performance(out1, measures = list(mlr::rmse, mlr::rae, mlr::rsq))

##         rmse          rae          rsq 
## 910.48403738   0.99480963  -0.01248682

Wyniki klasyfikacji na bodyfat:

performance(out2, measures = list(mlr::acc, mlr::auc, mlr::f1))

##       acc       auc        f1 
## 0.5700000 0.7524038 0.2950820

Jak możemy zobaczyć, wyniki dopasowania nie są powalające. Wręcz przeciwnie - przy klasyfikacji nasz model działa gorzej, niż gdyby zgadywał. W dalszej części będziemy starali się je poprawić. Najpierw jednak zobaczmy, jak wygląda maszyna dla zadania klasyfikacji na wykresie (patrzymy na obszary klasyfikacji w dwóch wymiarach - Density i Weight, dla reszty wymiarów biorąc wartości średnie ze zbioru treninigowego):

min(bodyfat[trainind2, c("Density")]) -> dmin
max(bodyfat[trainind2, c("Density")]) -> dmax
min(bodyfat[trainind2, c("Weight")]) -> wmin
max(bodyfat[trainind2, c("Weight")]) -> wmax
expand.grid(seq(dmin, dmax, length.out = 100), seq(wmin, wmax, length.out = 100)) -> denwei
colMeans(bodyfat[trainind2,-15]) -> means
matrix(rep(means,10000), ncol = 14, byrow = TRUE) -> rest
rest[,1] <- denwei[,1]
rest[,3] <- denwei[,2]
colnames(rest) <- colnames(bodyfat)[1:14]
predsgrid1 <- predict(model2, newdata = as.data.frame(rest))

ggplot(data = cbind(Response = predsgrid1$data$response, as.data.frame(rest)), aes(x=Density, y=Weight, color = Response)) +
  geom_point() +
  theme_bw() +
  geom_point(data = bodyfat[trainind2, ], aes(x=Density,y=Weight,shape=binaryClass), color = "black") +
  ggtitle("Klasyfikacja svm") +
  labs(shape = "Target")

Jak widzimy, maszyna klasyfikuje prawie wszystkie punkty danych jako “N”, natomiast obszar, w którym punkty oznaczałaby jako “P” jest bardzo wąski względem wymiaru Weight - wynika to z faktu, że wymiar ten posiada zupełnie inną skalę niż wymiar Density (różnica dwóch rzędów wielkości), podczas gdy wektory we wszyskich kierunkach mają podobne długości.

3 Skalowanie

Sprawdzimy teraz, jak poprawi się skuteczność modeli po przeskalowaniu danych. Stworzymy teraz dwa dodatkowe modele, w których przeskalujemy kolumny.

lrn1s <- makeLearner("regr.svm", par.vals = list(scale = TRUE))
lrn2s <- makeLearner("classif.svm", predict.type = "prob", par.vals = list(scale = TRUE))

model1s <- train(lrn1s, task1)
model2s <- train(lrn2s, task2)

out1s <- predict(model1s, newdata = apartments[testind1,])
out2s <- predict(model2s, newdata = bodyfat[testind2,])

cat("Wyniki regresji na apartments (ze skalowaniem): \n")

## Wyniki regresji na apartments (ze skalowaniem):

performance(out1s, measures = list(mlr::rmse, mlr::rae, mlr::rsq))

##        rmse         rae         rsq 
## 170.9862203   0.1898419   0.9642919

cat("\nWyniki klasyfikacji na bodyfat (ze skalowaniem): \n")

## 
## Wyniki klasyfikacji na bodyfat (ze skalowaniem):

performance(out2s, measures = list(mlr::acc, mlr::auc, mlr::f1))

##       acc       auc        f1 
## 0.9000000 0.9707532 0.8979592

Widzimy zdecydowaną poprawę - bez skalowania modele były użyteczne w niewielkim stopniu, natomiast po przeskalowaniu mają już bardzo dobrą skuteczność. Zobaczmy jeszcze wykres dla tych samych dwóch wymiarów, co poprzednio:

predsgrid2 <- predict(model2s, newdata = as.data.frame(rest))

ggplot(data = cbind(Response = predsgrid2$data$response, as.data.frame(rest)), aes(x=Density, y=Weight, color = Response)) +
  geom_point() +
  geom_point(data = bodyfat[trainind2, ], aes(x=Density,y=Weight,shape=binaryClass), color = "black") +
  theme_bw() +
  ggtitle("Klasyfikacja svm po przeskalowaniu") +
  labs(shape = "Target")

Wniosek: skalowanie svm jest przydatne.

4 Hiperparametry

Teraz spróbujemy dostroić hiperparametry do obu modeli. Skorzystamy do tego celu z metod optymalizacji udostępnianych przez pakiet mlrMBO. Pozwala on na utworzenie nowego typu obiektu, jaki możemy przekazać do funkcji tuneParams() z pakietu mlr - TuneControlMBO. Hiperparametry będą dobierane poprzez zewnętrzny model. Aby przyspieszyć obliczenia, będziemy korzystać z wielowątkowości.

lrn1t <- makeLearner("regr.svm", par.vals = c(list(scale = TRUE), parres1$x))
lrn2t <- makeLearner("classif.svm", predict.type = "prob", par.vals = c(list(scale = TRUE), parres2$x))

model1t <- train(lrn1t, task1)
model2t <- train(lrn2t, task2)

out1t <- predict(model1t, newdata = apartments[testind1,])
out2t <- predict(model2t, newdata = bodyfat[testind2,])

Wyniki crossvalidacji dla regresji:

parres1

## Tune result:
## Op. pars: kernel=radial; cost=6.05e+03; gamma=0.00217
## rmse.test.rmse=152.2930182

Wyniki na zbiorze testowym:

performance(out1t, measures = list(mlr::rmse, mlr::rae, mlr::rsq))

##        rmse         rae         rsq 
## 131.7179611   0.1540001   0.9776607

Wyniki crossvalidacji dla klasyfikacji:

parres2

## Tune result:
## Op. pars: kernel=linear; cost=0.201
## acc.test.mean=0.9605229

Wyniki na zbiorze testowym:

performance(out2t, measures = list(mlr::acc, mlr::auc, mlr::f1))

##       acc       auc        f1 
## 0.9600000 0.9951981 0.9591837

5 Analiza zależności od zmiennych

Na koniec przyjrzymy się wykresom, generowanym za pomocą pakietu DALEX, pozwalającym zrozumieć nam zachowanie modeli.

Na wszystkich wykresach “before” oznacza model przed dostrojeniem hiperparametrów, “after” - po dostrojeniu, a “compare” - model lasu losowego, dla porównania.

Najpierw wykresy dla regresji.

lrn1c <- makeLearner("regr.ranger")
model1c <- train(lrn1c, task1)
out1c <- predict(model1c, newdata = apartments[testind1,])

custom_predict1 <- function(object, newdata) {pred <- predict(object, newdata=newdata)
                                              response <- pred$data$response
                                              return(response)}
expl1bef <- explain(model1s, data=apartments[testind1,], 
                    y=apartments[testind1,"m2.price"], predict_function = custom_predict1, label="before")
expl1aft <- explain(model1t, data=apartments[testind1,], 
                    y=apartments[testind1,"m2.price"], predict_function = custom_predict1, label="after")
expl1com <- explain(model1c, data=apartments[testind1,], 
                    y=apartments[testind1,"m2.price"], predict_function = custom_predict1, label="compare")


mp1bef <- model_performance(expl1bef)
mp1aft <- model_performance(expl1aft)
mp1com <- model_performance(expl1com)

pdp1bef  <- variable_response(expl1bef, variable =  "construction.year", type = "pdp")
pdp1aft  <- variable_response(expl1aft, variable =  "construction.year", type = "pdp")
pdp1com  <- variable_response(expl1com, variable =  "construction.year", type = "pdp")

plot(mp1bef, mp1aft, mp1com)

plot(pdp1bef, pdp1aft, pdp1com)

Pierwszy wykres pokazuje rozkład wartości bezwzględnej błędów predykcji. Widać wyraźnie, że po dostrojeniu hiperparametrów reszty w pewnym stopniu zmalały.

Drugi wykres prezentuje reakcję modeli na poszczególne wartości zmiennej “construction.year”. Jak widać, model dostrojony reaguje silniej niż niedostrojony. Oba jednak prezentują zależność wielomianową, podczas gdy model drzewiasty jest raczej skokowy.

lrn2c <- makeLearner("classif.ranger", predict.type = "prob")
model2c <- train(lrn2c, task2)
out2c <- predict(model2c, newdata = bodyfat[testind2,])

custom_predict2 <- function(object, newdata) {pred <- predict(object, newdata=newdata)
                                              response <- pred$data[,3]
                                              return(response)}

y2 <- as.numeric(bodyfat[testind2,"binaryClass"]) -1

expl2bef <- explain(model2s, data=bodyfat[testind2,], 
                    y=y2, predict_function = custom_predict2, label="before")
expl2aft <- explain(model2t, data=bodyfat[testind2,], 
                    y=y2, predict_function = custom_predict2, label="after")
expl2com <- explain(model2c, data=bodyfat[testind2,], 
                    y=y2, predict_function = custom_predict2, label="compare")


mp2bef <- model_performance(expl2bef)
mp2aft <- model_performance(expl2aft)
mp2com <- model_performance(expl2com)

pdp2bef  <- variable_response(expl2bef, variable =  "Weight", type = "pdp")
pdp2aft  <- variable_response(expl2aft, variable =  "Weight", type = "pdp")
pdp2com  <- variable_response(expl2com, variable =  "Weight", type = "pdp")

plot(mp2bef, mp2aft, mp2com)

plot(pdp2bef, pdp2aft, pdp2com)

Ponownie, pierwszy wykres dla klasyfikatorów prezentuje rozkład reszt (błędem jest tu dla nas prawdopodobieństwo przeciwnej klasy). Ponownie, dostrojony model ma mniejsze błędy.

Drugi wykres pokazuje, że po dobraniu hiperparametrów, model przestał zwracać tak dużą uwagę na zmienną Weight, podobnie jak model drzewiasty.